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ABSTRACT 


The discontinuous Galerkin (DG) method has been accepted in the last decade 
by the geosciences community as an important component of geophysical fluid dy¬ 
namics. The high-order accuracy, geometric flexibility to use unstructured grids, 
local conservation, and monotonicity properties of the DG method make it a prime 
candidate for the construction of future ocean and shallow water models. 

This study focuses on formatting real bathymetry data of the Indian Ocean in 
order to simulate the propagation stage of the Indian Ocean tsunami that occurred 
on December 26, 2004, by using a DG model. In order to validate this simulation 
the study uses real measurements. The model results are compared to tide gauge 
data from several stations around the Indian Ocean, satellite altimetry, and held 
measurements. These results show that the model gives accurate estimates of arrival 
times in distant locations. 
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I. 


INTRODUCTION 


The 9.3 magnitude Sumatra-Andaman undersea earthquake occurred at 00:58:53 
UTC (07:58:53 local Indonesian time) on December 26, 2004. This earthquake gener¬ 
ated a megatsunami that was among the deadliest disasters in modern history, killing 
more than 226,000 people [1] and is known as the 2004 Indian Ocean tsunami. 

The modern science triggered by this and other similar disasters aimed to 
predict when and where an earthquake would occur and where a tsunami might be 
initiated. While the accurate prediction of an earthquake itself remains a major 
challenge, simulating the propagation and inundation of a tsunami wave is now a 
distinct possibility. Knowing where the faults are through the use of real bathymetry 
data, numerical simulations can be used to identify the regions that are likely to be 
affected by tsunami waves and establish seismic criteria for issuing tsunami warnings 
in the case of actual tsunamis. Furthermore, by using the existing technology of 
monitoring seismic activity, the simulation of the generated tsunami can give sufficient 
warning to the regions that are in the path of the wave. 

Tsunami simulation, as well as any other numerical modeling, requires the 
discretization of the physical space. The two major strategies used are regular and 
unstructured discretizations [2]. Models that are based on regular meshes and have 
been successfully applied to the 2004 Indian Ocean tsunami are as follows: the MOST 
(Method of Splitting Tsunami) propagation model by Titov and Synolakis [3, 4], 
which is based on a hnite-difference numerical approximation to the nonlinear shal¬ 
low water wave equations, and the fully nonlinear and dispersive Boussinesq model 
(FUNWAVE) by Watts [5]. An example of a model based on unstructured grids is 
the TsunAWI model by the tsunami modeling group of the Alfred-Wegener Institute 
[ 2 ]. 

The discontinuous Galerkin (DG) method, because of the high-order accuracy, 
geometric flexibility to use unstructured grids, local conservation, and monotonicity 
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properties, has been accepted in the last decade by the geosciences community as a 
critical component of the computational geophysical fluid dynamics and specihcally, 
in the construction of future ocean and shallow water models [6]. Schwanenberg 
and Kongeter [7] hrst used the DG method for the planar shallow water equations 
followed by the works of Li and Liu [8] as well as Ainzinger and Dawson [9]. Later, 
Dupont and Lin [10], Eskilsson and Sherwin [11], Remade et ah [12], and Kubatko 
et ah [13] constructed shallow water models on triangles using a collapsed local 
coordinate (modal) discontinuous Galerkin method. Giraldo et ah [14], hrst used the 
DG method for the shallow water equations on the sphere and later by construction 
on triangular domains [15]. Giraldo and Warburton et al. [6] also applied the method 
to the two-dimensional oceanic shallow water equations. Motivation for the choice 
to use unstructured grids is related to its ability to better represent the continental 
coastlines and simulate wave processes on different scales. A tsunami’s motion is 
characterized by relatively fast speed and low amplitude in the deep ocean, and when 
it reaches shallow waters the speed decreases and the amplitude increases. Near the 
shoreline, the tsunami steepens and may break. 

Giraldo and Warburton [6] DG formulation used high-order Langrange poly¬ 
nomials on the triangle using nodal sets up to the 15th order. By using six test cases 
(three of which had analytic solutions) they showed that the high-order triangular 
DG method exhibits exponential convergence. 

According to [16] the evolution of earthquake-generated tsunamis has three 
distinctive stages: generation, propagation, and run-up. This study, as being part 
of a research in regards to the development of a triangular discontinuous Galerkin 
oceanic shallow water model, focuses on formatting real bathymetry data of the In¬ 
dian Ocean in order to simulate the propagation stage of the Indian Ocean tsunami of 
December 26, 2004, by using this DG model. The initial conditions (i.e., the tsunami’s 
generation) are taken from the reconstruction of the rupture as described in [17] and 
based on [18], while the run-up stage is out of the scope for this thesis. Real measure- 
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ments are used for the validation of this simulation. The model results are compared 
to tide gauge data from several stations around the Indian Ocean, satellite altimetry 
(Jason-1, Topex-Poseidon and Envisat data), and held measurements. These results 
show that the model gives accurate estimates of arrival times in distant locations. 

This study is organized as follows. Chapter II describes the spatial discretiza¬ 
tion of the governing shallow water equations and the time integrator used in this 
model. Chapter III describes the real bathymetry data, the initial conditions used 
for the simulation, and the real measurement data used for comparisons. Chapter 
IV presents the results from the comparison and Chapter V discusses the conclusions 
and recommendations. Appendix A describes the derivation of the analytic solution 
for the Munk problem. This solution is given in order for the validation of the lateral 
diffusion operators in future studies. 
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II. 


BACKGROUND 


A. GOVERNING EQUATIONS 
1. Shallow Water Equations 

The oceanic shallow water equations are a system of nonlinear partial differen¬ 
tial equations which govern the motion of a viscous incompressible fluid in a shallow 
depth. The underlying assumption is that the depth of the fluid is small compared to 
the wave length of the disturbance. The Indian Ocean is not shallow since the depth 
is several kilometers, but the devastating tsunami on December 26, 2004, involved 
waves that were tens or hundreds of kilometers long, so the shallow water approxi¬ 
mation provides a reasonable model in this situation. A tsunami is the ideal shallow 
water flow problem because of the barotropic motion of the water column. 

The shallow water equations in conservation form are: 

^ + V . F(q) = S(q) (IM) 

where q = (0, are the conservation variables, 

F(q) ={ ] (IL2) 

y (l)u®u+ {(pu) j 

is the flux tensor, and 

^(g) = - \ (II.3) 

y / (fc X (j)u) + (j)'V(j)^ ~ ^ / 


is the source function where: 

- the nabla operator is dehned as V = dyY 

- ® denotes the tensor product operator 

- (p = gh is the geopotential height 

(where g is the gravitational constant and h is the free surface height of the fluid), 

- (p^ is the bathymetry 
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- u = {u, vY is the velocity vector 

■ f = fo + P{y ~ Um) is the Coriolis parameter 

- A; = (0, 0,1)^ is the unit normal vector of the x - y plane 

- the vector r is the wind stress 

- the constant 7 is the bottom friction 

- the term I 2 is the rank -2 identity matrix 

B. TRIANGULAR DISCONTINUOUS GALERKIN METHOD 

This section describes the discretization of the shallow water equations by the 
discontinuous Galerkin method on triangles, as dehned by Giraldo and Warburton 
[ 6 ], 

1. Basis Functions 

This method demands the decomposition of the domain hi into Ne non-overlapping 
triangular elements such that 

Ne 

e=l 

and the introduction of a nonsingular mapping a; = T(^) which dehnes a transforma¬ 
tion from the physical Cartesian coordinate system x = {x,yY to the local reference 
coordinate system ^ 7 )^, dehned on the reference triangle 

= {(C,h), -1 < < 1, ^ + V< 0,}. 

The local element-wise solution q can be represented by an iVth order poly¬ 
nomial in ^ as 

Mn 

<7Ar(«) = EV'.(«)<7K(«.) (11.4) 

i=l 

where represents M^r = -|- l)(iV -|- 2) interpolation points and are the 

associated multivariate Lagrange polynomials. 
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The Lagrange polynomials (nodal basis fnnctions), on the reference 

triangle are dehned as 




(2i + l){i + j + 1) ^0 0 


(W) 


1 ~ h V p2i+l,0 

2 J ^ 


iv) 


(11.5) 


where represents the nth order Jacobi polynomial in the interval —1 < ^ 1) 

k = i + j{N + 1) + 1, and the indices vary as 0 < i,j;i + j < N, and k = 1,..., M^. 
An explicit formnla for the Lagrange basis is 


Miv 

fc=l 

where the indices are dehned as i,j,k = 1,..., M^. By nsing the cardinal property 
of the Lagrange polynomials 

Mm 

^ij ^ ) ■^ik'^k {^ji hi) ) 
fc=l 

where S is the Kronecker delta fnnction, it can be written as 


Ak={^kHAvi)y 


(11.7) 


Recognizing that 


Vjk — <fk hi) 


( 11 . 8 ) 


is the generalized Vandermonde matrix and using Equations (IL6), (IL7), and (IL8), 
the Lagrange polynomials can be constructed as follows: 

j. 
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2. Integration 

For any two functions (/) and (g), the 2D (area) integration X 4 that is required 
by the weak formulation of Galerkin methods is given by 


. Me 

[/,g] = f{x)g{x)dx = ^wf I r(^i) I /(^i)g(^i) 

where Mq is a function of C which represents the order of the cubature approximation. 
For Wi and the high-order cubature rules for the triangle are used. 

The DG method requires the evaluation of boundary integrals. This is the 
mechanism by which the fluxes across element edges are evaluated and allows the dis¬ 
continuous elements to communicate. According to [ 6 ], the ID (boundary) integration 
Xg is given by 


^B[/,g] = j f{x)g{x)dx = I /(^i)g(^i) 

e i=0 

where Q represents the order of the quadrature approximation. The use of Q 
gives order 2N — 1 accuracy. 


= N 


3. Semi-discrete Equations 

The use of the discontinuous Galerkin discretization in Equation (II.1) gives 
the weak form of the DG method 

— i^TV ■ V — ilji{x)dx = ~ ^ n ■ F*^dx (II.10) 

where = F^q^q) and Sqq = S{qqq) with F and S given by Equations (II.2) and 
(II.3), respectively. 

In the boundary integral of Equation (II. 10) n is the outward pointing unit 
normal vector of the element edge Fg and F*^^ is the Rusanov numerical flux as given 
by the following equation: 



n = V 


2 L' 


Fn (qn) + Fn {q§) -\M{QN-qN 


n 


( 11 . 11 ) 


where A = max (^\U^\ + \ U^\ + with . n being the normal 

component of velocity with respect to the edge Fg, and the superscripts L and R 



represent the left and right sides of the element edge. The normal vector n is defined 
as pointing outward from left to right. 

Integrating Equation (II. 10) by parts once more yields the strong form 

Mx) + V-FM-SN^dx = n-{FN- F*j^) dx. ( 11 . 12 ) 

4. Matrix Form of the Semi-discrete Equations 

Using the polynomial approximation 

Mjv 

1=1 

the semi-discrete system can be written as 



ipiipj dx 


dt 


+ 



i/ji’Vi/jj dx ■ Fj 



Next, the elemental matrices M?- are dehned as follows: 


Fir 

( 11 , 13 ) 



which makes it possible to write Equation (11.13) in the matrix form: 

dq-j 


M. 




dt 


D\ 


ijj 


M‘Sj = (Mt,Y (F - F- 


( 11 . 14 ) 


where the superscript e denotes an element-based evaluation and s denotes a side- 
based (or edge-based) evaluation. 

The elemental matrices can be written, in terms of the computational variables 
(^), as follows: 

M!, = \r\ L = 

\J Q ig 



r 



di = \j%D%c + D2mi)i + \r\(Dl£ 


+ DlfQj 


M"- = 


r 


'ipi'ipjUdi = I + ntj 
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where the metric terms are dehned as |J®| = 2A® and |J^| = 2A^. The A® and A^ 
denote the area of element e and the length of edge s respectively. Note that the 
gradient operator in these integrals, V^, is dehned in the local reference coordinate 
system and hie and Tg denote the area and bonndary domains in the compntational 
space, namely the bonnds of integration for the master element. Thus, Equation 
(11.14) can be written in the following way: 


\r\M„^ + I J'l fj + {dUI + Dlnl) g‘i - 

( 11 , 15 ) 


= n‘jr-r), + nl(g‘-g‘), 


where + g^j. 

By dividing Equation (11.15) by | J®| and left multiplying by the following 
is obtained 


dqt 

dt 


+ [Die +Ate) fi + (Sly;+Ati: 




st = 


\J-\ 




(11.16) 


where the matrices are dehned as 


where 


D% = Dl = Mg'Dl, Mt, = Mg^Ml 

Me Mq 

^ij = Wk'tpiki’jk, Wki’ik'tpjk, 

k=l k=l 


( 11 , 17 ) 


( 11 , 18 ) 


r,« , aV'jt ™ 

Dij = 2^Wk'ipik-^, D^j = }_^WkiJik-^- (11.19) 

k=i k=i 

The Me and Mq denote the number of cubature (two-dimensional) and quadrature 
(one-dimensional) integration points required to achieve order 2N — 1 accuracy, and 
ipik represents the function tjj at the i = 1, interpolation points evaluated at 

the integration point k. 
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Since the mass matrix is constant (i.e., not a function of x) then the use of 
Equations (11.18) and (11.19) gives the following: 

^ij = Wk^ik'ipjk, Dij = ^ WkAk^!^, 

k=l k=l k=l 

where 

A = 

is the basis function premultiplied by the inverse mass matrix. Absorbing the mass 
matrix within the test function ijj completely eliminates the mass matrix from Equa¬ 
tion (11.16) without making any approximations. 


5. Time Integrator 

In order to advance the solution in time while retaining high-order accuracy, 
the strong stability preserving (SSP) Runge-Kutta third-order RK3 method (see [6]) 
is used. The semi-discrete (in space) equations are written as follows: 


dq 


= S{q). 


The SSP temporal discretization of this vector equation is 


fork = 1,..., 3 : 

q^ = a’^q^+ a’lq’^-^ + p’^AtSiq^-^) 


where q^ = g"", q^ = and the coefficients a and f3 are given in Table I. 



k=i 

k=2 

k=3 

ao 

1 

3/4 

1/3 

(y.\ 

0 

1/4 

2/?> 


1 

1/4 

2/3 


Table I. Coefficients for the Strong Stability Preserving - Third Order Runge-Kutta 
Method. 
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III. SIMULATION OF THE INDIAN OCEAN 

TSUNAMI 

A. MODEL DESCRIPTION 

1. Boundary Conditions 

As was mentioned above, this thesis study focuses on the simulation of the 
propagation stage of the Indian Ocean tsunami while the inundation modeling is out 
of its scope. For this reason, the shoreline is being treated as a hxed-wall boundary. 
These no-flux boundary conditions are enforced at all land boundaries by virtue of 
the statement 

n ■ = 0 (III-l) 

which eliminates the normal component of the velocity to the no-flux boundary with¬ 
out altering the tangential component (also known as free slip boundary conditions). 
Furthermore, the points near the coastlines with depth less than 11 meters are con¬ 
sidered to have 11 meters depth. 

Since the model does not currently have inundation algorithms, it only makes 
sense to compare observations with the model results using arrival times and heights 
of the first waves. This way the model results will not be affected by the errors 
associated with representation of boundary conditions along the coastlines. 

2. Tides / Coriolis / Viscosity 

This study is based on two experiments. Each experiment had ten hours of 
simulations with time steps of 1.5 seconds. During the simulations, tidal forcing 
and the Coriolis effect were not included. In the hrst experiment, the value of the 
horizontal viscosity coefficient was v = 1000, while in the second the value was zero. 
Both simulations yielded the same results. This can be explained by using scale 
analysis for the shallow water momentum equations, which indicates that the relation 
of the horizontal viscosity to the rest of the forces is 10“®, on the order of machine 
single precision. 
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B. REAL BATHYMETRY DATA 


For this study, real bathymetry data of the Indian Ocean was provided by the 
Alfred-Wegener Institute for Polar and Marine Research (tsunami modeling group). 
The unstructured mesh was produced with the mesh generator TRIANGLE by Jonathan 
Shewchuk [20]. As described in [2], the resulting meshes were smoothed to improve 
the overall mesh quality. The coarsest resolution of the mesh is 14 kilometers in the 
deep ocean, 500 meters in the coastal areas and in areas with high bathymetry gra¬ 
dients (canyons), and less than 100 meters in the northern tip of Sumatra. The total 
number of elements used was 130,445 created by 66, 715 grid points. Figure 1 presents 
the grid points used for this simulation study and the Indian Ocean bathymetry in 
meters. 



30°E 45°E 60°E 75°E 90°E lOS^E 

Longitude 



Longitude 


4000 

3000 

2000 

1000 

0 

-1000 

-2000 

-3000 

-4000 

-5000 


Figure 1. Mesh Grid Points Provided by AWI and Indian Ocean Bathymetry 


The conversion of the above real bathymetry data from latitude/longitude co¬ 
ordinates into Gartesian coordinates (x,y) was made by using the MATLAB function 
m-ll2xy as described in [21]. The used projection was a Mercator projection with 
the properties given in Table II. 
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Min Longitude 

Max Longitude 

Min Latitude 

Max Latitude 

030°U 

110 °U 

35°^ 

35°N 


Table II. Properties of the Used Mercator Projection 


C. TSUNAMI SOURCE MODEL 

The main cause of tsunami waves is the abruptly motion of converging or 
destructive plate boundaries which vertically displace the overlying water. The Indian 
Ocean tsunami was generated by the static sea floor uplift caused by an abrupt slip 
of the India/Burma plate interface. This permanent, vertical sea floor displacement 
can be computed using the static dislocation formulae from [18]. Inputs for this 
computation are the fault plane location, its depth, strike, dip, slip, length, width, 
seismic moment, and rigidity. The initial conditions for the tsunami simulation were 
taken from the reconstruction of the rupture as described in [17]. According to this 
study, the fault extent constrained by observed tsunami arrival time to the northwest, 
east and south of the slip zone indicates a fault zone of approximately 1000 kilometers 
by 200 kilometers. The epicenter location lies on the southern end of the fault zone. 
To accommodate trench curvature, this fault plane was broken into two segments as 
depicted in Figure 2. The points that form these planes and used for this study are 
presented in Table III. 
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Figure 2. Initial Conditions for the Simulation of the December 26,2004, Indian Ocean 
tsunami [After: [17]] 


Fault Segment 

Corner 

Latitude 

Longitude 

Sea surface 
elevation 

Northern Fault Segment 

NW 

11.78°A^ 

092.12°E 

+5.07 

m 

Northern Fault Segment 

sw 

5.6°N 

093.22°E 

+5.07 

m 

Northern Fault Segment 

NE 


094.02°E 

-4.75 

m 

Northern Fault Segment 

SE 

6 .0° AT 

094.97°E 

-4.75 

m 

Southern Fault Segment 

NW 

5.33°Ar 

093.25°E 

+5.07 

m 

Southern Fault Segment 

SW 

2.97° AT 

094.35°E 

+5.07 

m 

Southern Fault Segment 

NE 

6 .0° AT 

094.97°E 

-4.75 

m 

Southern Fault Segment 

SE 

3.88° A^ 

095.97°E 

-4.75 

m 


Table III. Points That Form the Two Fault Segments [After: [17]] 
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D. REAL MEASUREMENT DATA 
1. Tide Gauge Records 

Tide gauge records are often used to provide valuable information on the 
tsunami arrival times and the changes in sea level due to tsunamis. After re¬ 
moving the tidal effect these data provide the changes in water level due to tsunamis 
alone. 

The Indian Ocean tsunami was recorded on all the tide gauges located in the 
Indian Ocean. A comprehensive overview of these stations and an analysis of the 
tsunami records is published in [22]. For this thesis study tide gauge records from 
hfteen stations in the Indian Ocean were used. The locations of these stations are 
shown in Figure 3. In order to increase the accuracy of the results, each station’s 
location was interpolated using the three grid points of the triangular element in 
which the station belongs. The tsunami characteristics of each station’s records are 
shown in Table IV. 
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Figure 3. Locations of the Tide Gauges Where the December 26, 2004, Tsunami 
Waves Were Recorded [After: [22, 23, 25]] 
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No 

Station 

Coordi¬ 

nates 

Sampling 

interval 

(min) 

First 

Arrival 
time (UTC) 


wave 

Travel 

time 

1 

Paradip 

20.26°V 

086.70°E 

6 

03 : 27 

2 

hrs 28 min 

2 

Vishakhapatham 

17.65°V 

083.28°E 

5 

03 : 35 

2 

hrs 36 min 

3 

Chennai 

13.10°V 

080.32°E 

5 

03 : 33 

2 

hrs 34 min 

4 

Tuticorin 

08.75°V 

078.20°E 

6 

04 ; 24 

3 

hrs 25 min 

5 

Kochi (Cochin) 

09.97°V 

076.27°E 

6 

05 : 41 

4 

hrs 42 min 

6 

Mormugao 

15.42°V 

073.80°E 

5 

06 : 53 

5 

hrs 54 min 

7 

Okha 

22.50°V 

069.10°E 

6 

09 : 03 

8 

hrs 04 min 

8 

Colombo 

06.93°V 

079.83°E 

2 

03 : 49 

2 

hrs 50 min 

9 

Hanimaadhoo 

06.77°V 

073.18°E 

2 

04 : 30 

3 

hrs 31 min 

10 

Male 

04.18°V 

073.52°E 

4 

04 : 14 

3 

hrs 15 min 

11 

Gan 

00 .68°A 

073.17°E 

4 

04 ; 16 

3 

hrs 17 min 

12 

Diego Garcia 

07.30°A 

072.38°E 

6 

04 : 45 

3 

hrs 46 min 

13 

Port Louis 

20.15°^ 

057.50°E 

2 

07 : 46 

6 

hrs 47 min 

14 

Salalah 

17.00°V 

054.00°E 

4 

08 : 08 

7 

hrs 09 min 

15 

Pointe La Rue 

04.68°A 

055.53°E 

4 

08 : 16 

7 

hrs 17 min 


Table IV. Tsunami Characteristics Estimated From the Indian Tide Gauge Records 
[After: [22]] 
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The records of the hrst seven stations of Table IV (these stations are depicted 
with their names in black color in Figure 3) were obtained from the National Institute 
of Oceanography (NIO), Goa, India [23]. These tide gauges are part of the network 
maintained by the Survey of India (SOI) agency along the coast of India. According 
to [24], SOI tide gauges were either mechanical float-type analog gauges or pressure¬ 
sensor gauges. The analog records from Vishakhapatham, Chennai, and Mormugao 
were digitized by SOI at an interval of 5 minutes, while those from Okha at an interval 
of 6 minutes. The pressure-sensor gauges from Paradip, Tuticorin, and Kochi had 
sampling intervals of 6 minutes. All these are high-quality, de-tided records. Some 
extensive gaps in the data were due to instruments and transmission problems. A 
typical example of these records is presented in Figure 4. 


Paradip - Observation Data - Pos.: 20.26N ; 086.70E 



Figure 4. Time Series of the Sea Surface Elevation in Paradip - India [After: [23]] 
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The records of the next eight stations of Table IV (these stations are depicted 
with their names in red color in Figure 3) were obtained from the University of 
Hawaii’s Sea Level Center database (Honolulu) [25]. These are digital Global Sea 
Level Observing System (GLOSS) stations with 2 minutes sampling intervals for 
Colombo, Hanimaadhoo, and Port Louis; 4 minutes for Male, Gan, Salalah, and 
Pointe La Rue; and 6 minutes for Diego Garcia. Most of these stations are located on 
isolated islands and therefore the records were not signihcantly affected by reflections 
[22]. These records were real-time sea level measurements (that included tides) and 
were de-tided by using the t^tide program as described in [29] and [30]. Two problems 
during this procedure were the extensive gaps in the data due to instruments and 
transmission problems and the limited total record time (only one week). The results 
are presented in Figures 5 through 12. In all the above cases, except for Diego Garcia, 
the tidal effect was signihcantly reduced. 


Colombo - Observation Data - Pos.: 06.93N ; 079.83E Colombo - De-Tidal Data - Pos.; 06.93N ; 079.83E 




Figure 5. Results From De-tidal Procedure for Colombo - Sri Lanka Tide Gauge 
Records. [After: [25],[29],[30]] 
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200 


Hanimaadhoo - Observation Data - Pos.: 06.77N ; 073.18E 



Date 


Hanimaadhoo - De-Tidal Data - Pos.; 06.77N ; 073.18E 



Figure 6. Results From De-tidal Procedure for Hanimaadhoo - Maldives Tide Gauge 
Records. [After: [25],[29],[30]] 


Male - Observation Data - Pos.; 04.18N ; 073.52E 



Male (De-Tidal) - Pos.: 04.18N ; 073.52E 



Figure 7. Results From De-tidal Procedure for Male - Maldives Tide Gauge Records. 
[After: [25],[29],[30]] 
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Gan - Observation Data - Pos.: 00.68S ; 073.17E 



Gan - De-Tidal - Pos.: 00.68S ; 073.17E 


80 
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-40 

-60 


.801 - 1 -^^- 1 -^^- 1 
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Figure 8. Results From De-tidal Procedure for Gan - Maldives Tide Gauge Records. 
[After: [25],[29],[30]] 


Diego Garcia - Observation Data - Pos.: 07.30S : 072.38E Diego Garcia - De-Tidal Data - Pos.: 07.30S ; 072.38E 




Figure 9. Results From De-tidal Procedure for Diego Garcia - UK Tide Gauge 
Records. [After: [25],[29],[30]] 
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Port Louis - Observation Data - Pos.: 20.15S ; 057.50E 



Port Louis - De-Tidal Data - Pos.: 20.15S ; 057.50E 





25 26 27 28 29 30 31 01 

Date 


Figure 10. Results From De-tidal Procedure for Port Louis - Mauritius Tide Gauge 
Records. [After: [25],[29],[30]] 


Salalah - Observation Data - Pos.: 17.00N ; 054.00E Salalah - De-Tidal Data - Pos.: 17.00N ; 054.00E 




Figure 11. Results From De-tidal Procedure for Salalah - Oman Tide Gauge Records. 
[After: [25],[29],[30]] 
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Pointe La Rue - Observation Data - Pos.: 04.68S ; 055.53E Pointe La Rue - De-Tidal Data - Pos.; 04.68S ; 055.53E 




Figure 12. Results From De-tidal Procedure for Pointe La Rue - Seycelles Tide Gauge 
Records. [After: [25],[29],[30]] 


2. Satellite Altimetry 

Radar altimeters on board the Jason-1, TOPEX/Poseidon, and Envisat satel¬ 
lites obtained profiles of sea surface height on trajectories across the Indian Ocean 
between two and four hours after the Sumatra earthquake on December 26, 2004. 

Jason-1 crossed the equator at 2:55 UTC, approximately two hours after the 
earthquake [26]. The data were received hours to days after ’’real time” so it was 
too late to detect and warn about the upcoming tsunami, but it can be used to 
validate new models. The data used for this study were obtained from NOAA [27] 
and originated from the satellite crossings of Figure 13 and Table V. 
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Figure 13. Satellites Crosses During the 2004 Indian Ocean Tsunami 


No 

Satellite 

Time after the Earthquake 

Color in Figure 13 

1 

Jasonl 

02 hrs 00 min 

Red 

2 

TOPEX / Poseidon 

02 hrs 05 min 

Black 

3 

Envisat 

03 hrs 15 min 

Green 


Table V. Profiles of Sea Surface Height Obtained by Satellites During the 2004 Indian 
Ocean Tsunami [After: [27]] 
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3. Field Measurements 

During the Indian Ocean tsunami a Belgian yacht, the Mercator, was anchored 
about 1.6 kilometers off the Phuket coast (Thailand) at 07.715°iV, 098.28°£' (Figure 
14). The yacht’s depth gauge was operational and measured changing wave heights 
during the tsunami [28]. Based on these measurements, the three main tsunami waves 
had trough-to-crest wave heights of 6.6, 2.2 and 5.5 meters. The first wave (trough) 
struck the yacht’s location at 02:38 UTC (1 hour and 39 minutes). 


Mercator - Pos.: 07.715S 098.28E 



Figure 14. The Mercator’s Anchored Location 
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IV. 


SIMULATIONS RESULTS 


A. COMPARISON TO TIDE GAUGE RECORDS 

The simulation results are separated into the following two groups. The hrst 
group represents simulations of sea surface height measurements at the tropical Indian 
Ocean (Male, Gan, Diego Garcia, Port Louis, and Pointe La Rue), Oman (Salalah) 
and west India tide gauge stations (Okha and Mormugao). The second group rep¬ 
resents simulations of sea surface height measurements at the east Indian tide gauge 
stations (Paradip, Vishakhapatham, Ghennai, Tuticorin, Kochi, and Golombo) as 
well as at Hanimaadhoo (Maldives). 

1. Tide Gauge Stations Located at the Tropical Indian 
Ocean and West of India Regions 

All these stations, because of their position, were directly affected by the 
tsunami waves emanating from the source area. 

a. Male Station (Maldives) 

This station was located at 04.18°iV, 073.52°i? and was approximately 
1,180 nautical miles from the west side of the north fault segment (Figure 15). The 
GLOSS tide gauge of this station was sampled in 4 minutes intervals. 

Figure 16 illustrates the agreement between the tide gauge record and 
the model results about the arrival time (3 hours and 15 minutes after the earthquake) 
and the sign (positive - wave crest) of the hrst wave. There is a difference (100 instead 
of 125 centimeters) on its maximum amplitude. 
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Latitude 


Male - Pos.: 04.18N 073.52E 



Figure 15. Male Tide Gauge Position 
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Sea surface elevation [cm] 


Arrival Time 


Tide gauge record Male vs. Model 



Figure 16. Comparison of arrival times - Male Tide Gauge vs. Model 
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b. Gan Station (Maldives) 

This station was located at OO.GS®^', 073.17°i? and was approximately 


1,260 nautical miles from the west side of the south fault segment (Figure 17). The 
GLOSS tide gauge of this station was sampled in 4 minutes intervals. 

From Figure 18 it is obvious that the model correctly predicts the ar¬ 
rival time (3 hours and 17 minutes after the earthquake), the sign (positive - wave 
crest), and almost the amplitude of the hrst wave (60 instead of 75 centimeters). 

c. Diego Garcia Station (UK) 

This station was located at 07.30°S', 072.38°i7 and was approximately 
1,460 nautical miles from the west side of the south fault segment (Figure 19). The 
GLOSS tide gauge of this station was sampled in 6 minutes intervals. 

Figure 20 documents the agreement between the tide gauge record and 
the model results about the arrival time (3 hours and 46 minutes after the earth¬ 
quake) and the sign (positive - wave crest) of the hrst wave, although there is a small 
difference (35 instead of 42 centimeters) on its maximum amplitude. 

d. Port Louis Station (Mauritius) 

This station was located at 20.15°S', 057.50°i7 and was approximately 
2,600 nautical miles from the west side of the south fault segment (Figure 21). The 
GLOSS tide gauge of this station was sampled in 2 minutes intervals, ceased to operate 
for 1 hour [22], and the sign of the hrst wave was unknown. 

Figure 22 depicts the agreement between the tide gauge record and the 
model on the arrival time (6 hours and 47 minutes after the earthquake) of the hrst 
wave. According to the model, the hrst wave was positive (wave crest). 
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Latitude 


Gan - Pos.: 00.68S 073.17E 



Figure 17. Gan Tide Gauge Position 
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Arrival Time 

Tide gauge record Gan vs. Model 



Figure 18. Comparison of arrival times - Gan Tide Gauge vs. Model 


34 



































Latitude 


Diego Garcia - Pos.: 07.SOS 072.38E 



Figure 19. Diego Garcia Tide Gauge Position 


35 

















Sea surface elevation [cm] 


Arrival Time 

Tide gauge record Diego Garcia vs. Model 



Figure 20. Comparison of arrival times - Diego Garcia Tide Gauge vs. Model 
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Latitude 


Port Louis - Pos.: 20.15S 057.50E 



Figure 21. Port Louis Tide Gauge Position 
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Arrival Time 

Tide gauge record Port Louis vs. Model 



Figure 22. Comparison of arrival times - Port Louis Tide Gauge vs. Model 
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e. Salalah Station (Oman) 

This station was located at 17.00°iV, 054.00°i7 and was approximately 
2,550 nautical miles from the west side of the north fault segment (Figure 23). The 
GLOSS tide gauge of this station was sampled in 4 minutes intervals. 

Figure 24 depicts the general agreement between the tide gauge records 
and the model. 


/. Points La Rue Station (Seychelles) 

This station was located at 04.68°S', 055.53°i7 and was approximately 
2,350 nautical miles from the west side of the south fault segment (Figure 25). The 
GLOSS tide gauge of this station was sampled in 4 minutes intervals. 

From Figure 26 it is obvious that the model correctly predicts the ar¬ 
rival time (7 hours and 04 minutes after the earthquake), the sign (positive - wave 
crest), and nearly the amplitude of the hrst wave (50 instead of 80 centimeters). 

g. Mormugao Station (India) 

This station was located at 15.42'^iV —073.80*^77 and was approximately 
1,300 nautical miles from the west side of the north fault segment (Figure 27). The 
GLOSS tide gauge of this station was sampled in 5 minutes intervals. 

Figure 28 depicts the agreement on the arrival time (5 hours and 54 min 
after the earthquake) and the sign (positive - wave crest) of the hrst wave between the 
tide gauge record and the model. Also, the maximum amplitude of the hrst arrival 
wave (almost 55 centimeters) was accurately predicted. 
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Latitude 


Salalah - Pos.: 17.00N 054.00E 



Figure 23. Salalah Tide Gauge Position 
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Sea surface elevation [cm] 


Arrival Time 


Tide gauge record Salalah vs. Model 



Figure 24. Comparison of arrival times - Salalah Tide Gauge vs. Model 
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Latitude 


Pointe La Rue - Pos.: 04.68S 055.53E 



Figure 25. Pointe La Rue Tide Gauge Position 
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Sea surface elevation [cm] 


Arrival Time 

Tide gauge record Pointe La Rue vs. Model 



Figure 26. Comparison of arrival times - Pointe La Rue Tide Gauge vs. Model 
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Latitude 


Mormugao - Pos.: 15.42N 073.80E 



Figure 27. Mormugao Tide Gauge Position 
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Sea surface elevation [cm] 


Arrival Time 

Tide gauge record Mormugao vs. Model 



Figure 28. Comparison of arrival times - Mormugao Tide Gauge vs. Model 
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h. Okha Station (India) 

This station was located at 22.50°iV, 069.10°i? and was approximately 
1,900 nautical miles from the west side of the north fault segment (Figure 29). The 
GLOSS tide gauge of this station was sampled in 6 minutes intervals. 

The negative values of the observations in Figure 30 can be attributed 
to the tidal effect, which was not completely removed (see Figure 31). Taking into 
account this dehciency of data, the model correctly predicts the sign (positive - wave 
crest), arrival time (8 hours and 04 minutes) and amplitude of the hrst wave (8 cen¬ 
timeters). 
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Latitude 


Okha-Pos.:22.50N 069.10E 



Figure 29. Okha Tide Gauge Position 
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Sea surface elevation [cm] 


Arrival Time 

Tide gauge record Okha vs. Model 



Figure 30. Comparison of arrival times - Okha Tide Gauge vs. Model 
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Sea surface elevation [cm] 


Okha - Observation Data - Pos.: 22.SON ; 069.1 OE 



Figure 31. Time Series of the Sea Surface Elevation in Okha - India [After: [23]] 
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A summary of the simulation results for the tide gauge stations located 


at tropical Indian Ocean and west of India regions are presented in Table VI. 


No 

Tide Gauge 
Station/ 

Field measurement 

Coordinates 

Arrival Time 
Difference 
(min) 

Max wave’s height 
Relative error 

1 

Male 

04.18°iV 

073.52°F 

0 

-0.2 

2 

Gan 

00.68°^ 

073.17°F 

0 

-0.2 

3 

Diego Garcia 

07.30°A 

072.38°E 

0 

-0.16 

4 

Port Louis 

20.15°^ 

057.50°E 

0 

(?) 

5 

Salalah 

17.00°iV 

054.00°E 

0 

-0.1 

6 

Pointe La Rue 

04.68°A 

055.53°.F 

0 

-0.35 

7 

Mormugao 

15.42^iV 

073.80°E 

0 

-0.09 

8 

Okha 

22.50°iV 

069.10°E 

0 

0 


Table VI. Summary of Results for Tide Gauge Stations Located at the Tropical Indian 
Ocean and west of India regions 


Comments: 

(1) The tide gauge at Port Louis ceased to operate for one hour [22] 

(2) The relative error is evaluated as follows: Error = 

^ ^ hobs 
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2. Tide Gauge Stations Located South-southeast of In¬ 
dia 

The tsunami waves that arrived at the region of these stations emanating from 
the source (traveled in a parallel direction to the earthquake rupture) and to the in¬ 
fluence of tsunami waves reflected from the southwestern coasts of the Indonesian 
Islands. 


a. Paradip Station (India) 

This station was located at 20.26°iV, 086.70°i? and was approximately 
600 nautical miles from the north side of the north fault segment (Figure 32). The 
GLOSS tide gauge of this station was sampled in 6 minutes intervals. 

Figure 33 illustrates the agreement with respect to the arrival time 
(2 hours and 28 minutes after the earthquake) and the positive sign (wave crest) 
of the hrst wave between the tide gauge record and the model. However, there is 
a signihcant difference (25 instead of 125 centimeters) in the maximum amplitudes 
between the observed height and those predicted by the model. The waves that arrived 
at this region represent a combination of waves emanated from the source (traveling 
parallel to the earthquake rupture) and waves reflected from the southwestern coasts 
of the Indonesian Islands. Figures 34 and 35 are snapshots of the model’s sea surface 
elevation two and three hours after the earthquake, respectively. These hgures depict 
the waves generated by reflections from the southwestern coasts of the Indonesian 
Islands. Since this model does not have inundation algorithms, it cannot be expected 
to accurately reproduce the maximum height of the tide gauge records. 
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Latitude 


Paradip - Pos.: 20.26N 086.70E 



Figure 32. Paradip Tide Gauge Position 
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Sea surface elevation [cm] 


Arrival Time 

Tide gauge record Paradip vs. Model 



Figure 33. Comparison of arrival times - Paradip Tide Gauge vs. Model 
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02 hrs 00 min 



Longitude 


Figure 34. Model’s Sea Surface Elevation Two Hours After the Earthquake 
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Figure 35. Model’s Sea Surface Elevation Three Hours After the Earthquake 
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b. Vishakhapatham Station (India) 

This station was located at 17.65°iV, 083.28°i7 and was approximately 


700 nautical miles from the north side of the north fault segment (Figure 36). The 
GLOSS tide gauge of this station was sampled in 5 minutes intervals. 

From Figure 37 it is obvious that the model correctly predicts the pos¬ 
itive sign (wave crest) of the hrst wave but an earlier arrival time (2 hours and 26 
minutes instead of 2 hours and 36 minutes after the earthquake), and underestimates 
the amplitude of the first wave (75 instead of 125 centimeters). As in Paradip case, 
this region was affected by a combination of waves from the source and waves reflected 
from the southwestern coasts of the Indonesian Islands. 

c. Chennai Station (India) 

This station was located at 13.10°A^, 080.32°i7 and was approximately 
750 nautical miles from the west side of the north fault segment (Figure 38). The 
GLOSS tide gauge of this station was sampled in 5 minutes intervals. 

Figure 39 illustrates the agreement on the arrival time (2 hours and 
34 minutes after the earthquake) and the positive sign (wave crest) of the hrst wave 
between the tide gauge record and the model. The model, however, overestimates the 
maximum amplitude (125 instead of 70 centimeters). 
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Latitude 


Vishakhapatham - Pos.: 17.65N 083.28E 



Figure 36. Vishakhapatham Tide Gauge Position 
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Arrival Time 


Tide gauge record Vishakhapatham vs. Model 



Figure 37. Comparison of arrival times - Vishakhapatham Tide Gauge vs. Model 



















































Latitude 


Chennai - Pos.: 13.10N 080.32E 



Figure 38. Chennai Tide Gauge Position 
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Sea surface elevation [cm] 


Arrival Time 

Tide gauge record Chennai vs. Model 



Figure 39. Comparison of arrival times - Chennai Tide Gauge vs. Model 
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d. Tuticorin Station (India) 

This station was located at 08.75°iV, 078.20°i? and was approximately 


900 nautical miles from the west side of the north fault segment Figure(40). The 
GLOSS tide gauge of this station was sampled in 6 minutes intervals. 

From Figure 41 it is obvious that the model correctly predicts the pos¬ 
itive sign (wave crest) of the hrst wave but an earlier arrival time (3 hours and 15 
minutes instead of 3 hours and 25 minutes after the earthquake), and underestimates 
the amplitude of the hrst wave (50 instead of 90 centimeters). The earlier arrival time 
of this case can be explained by the fact that the arrival waves, during their path 
from the initialization to the position of the gauge, had to pass through areas with 
depths less than 11 meters. Due to the boundary conditions used in the simulations 
(all positions with depth less than 11 meters were treated as having a depth of 11 
meters in order to avoid negative water depth, this is due to the absence of inundation 
algorithms in the current version of the model) the model predicts earlier arrival times. 

e. Kochi Station (India) 

This station was located at 09.97°iV, 076.27°i7 and was approximately 
1,050 nautical miles from the west side of the north fault segment (Figure 42). The 
GLOSS tide gauge of this station was sampled in 6 minutes intervals. 

Figure 43 shows that the model correctly predicts the positive sign 
(wave crest) of the hrst wave but an earlier arrival time (4 hours and 22 minutes 
instead of 4 hours and 42 min after the earthquake), and overestimates the amplitude 
of the hrst wave (85 instead of 70 centimeters). The earlier arrival time can be 
explained again because of the no-hux boundary conditions used along the coastlines. 
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Latitude 


Tuticorin - Pos.: 08.75N 078.20E 



Figure 40. Tuticorin Tide Gauge Position 
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Tide gauge record Tuticorin vs. Model 



Figure 41. Comparison of arrival times - Tuticorin Tide Gauge vs. Model 
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Latitude 


Kochi - Pos.: 09.97N 076.27E 



Figure 42. Kochi Tide Gauge Position 
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Tide gauge record Kochi vs. Model 



Figure 43. Comparison of arrival times - Kochi Tide Gauge vs. Model 
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/. Colombo Station (Sri Lanka) 

This station was located at 06.93°iV, 079.83°i? and was approximately 
800 nautical miles from the west side of the north fault segment (Figure 44). The 
GLOSS tide gauge of this station was sampled in 2 minutes intervals, was damaged 
by the hrst tsunami wave, and did not operate for 5 hours and 40 minutes [22], 
According to eyewitness accounts, the second wave was the highest. 

Figure 45 depicts the agreement between the tide gauge record and the 
model that the hrst wave was positive (wave crest) and the highest wave was the 
second. There is a difference on the arrival time (2 hours and 35 min instead of 2 
hours and 50 min after the earthquake) for the hrst wave and also a diherence in its 
maximum amplitude (80 instead of 200 centimeters). 

g. Hanimaadhoo Station (Maldives) 

This station was located at 06.77°iV, 073.18°i7 and was approximately 
1,200 nautical miles from the west side of the north fault segment (Figure 46). The 
GLOSS tide gauge of this station was sampled in 2 minutes intervals. 

Figure 47 depicts the agreement between the tide gauge record and the 
model on that the hrst wave was positive (wave crest). There is a small diherence on 
arrival time (3 hours and 26 minutes instead of 3 hours 31 minutes after the earth¬ 
quake) for the hrst wave and also a diherence in its maximum amplitude (60 instead 
of 160 centimeters). 
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Latitude 


Colombo - Pos.: 06.93N 079.83E 



Figure 44. Colombo Tide Gauge Position 
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Tide gauge record Colombo vs. Model 



Figure 45. Comparison of arrival times - Colombo Tide Gauge vs. Model 
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Latitude 


Hanimaadhoo - Pos.: 06.77N 073.18E 



Figure 46. Hanimaadhoo Tide Gauge Position 


69 
















Sea surface elevation [cm] 


Arrival Time 

Tide gauge record Hanimaadhoo vs. Model 



Figure 47. Comparison of Arrival Times Hanimaadhoo Tide Gauge vs. Model 
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A summary of simulation’s results for the tide gauge stations located 


at South - Southeast of India is presented in Table VII. 


No 

Tide Gauge 
Station/ 

Field measurement 

Coordinates 

Arrival Time 
Difference 
(min) 

Max wave’s height 
Relative error 

1 

Paradip 

20.26°iV 

086.70°F 

0 

-0.8 

2 

Vishakhapatham 

17.65°iV 

083.28°F 

-10 

-0.4 

3 

Chennai 

13.10°iV 

080.32°F 

0 

0.78 

4 

Tuticorin 

08.750iV 

078.20°E 

-10 

-0.44 

5 

Kochi (Cochin) 

09.97°iV 

076.27°E 

-20 

0.38 

6 

Colombo 

06.93°iV 

079.83°E 

-15 

-0.6 

7 

Hanimaadhoo 

06.77°iV 

073.18°.F 

-5 

-0.6 


Table VII. Summary of Results for the Tide Gauge Stations Located at the South- 
Southeast of India regions 


Comment: 

(1) The relative error is evaluated as follows: Error- = ’^-^odei-Kbs 

^ ' hobs 
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B. COMPARISON TO SATELLITE ALTIMETRY 

1. Jason 1 

One of the three satellites used in this study was Jason—1. This satellite 
crossed the equator at 2:55 a.m. UTC, which was two hours after the earthquake 
[26], and its altimeter recorded the differences in the sea surface elevation in the 
Bay of Bengal. Figure 48 presents the model’s sea surface elevation two hours after 
the earthquake. The black line in this figure represents the points used to compare 
sea surface elevation in Figure 49 between the satellite’s altimetry data (blue line) 
and the model’s data (red dots). The model estimates very accurately the leading 
wave crest at about 5°S and the double peak structure between this and the equator, 
although the heights are lower than the observations. However, between 5°N to 12‘’N 
the model’s crest is in contrast to the observation’s trough. Finally, between 12‘^N to 
20°N the model and satellite data are in total agreement and specifically, the model 
estimates very accurately the leading wave crest at about 20°N. 


Sea surface elevation 2 hours after the earthquake 





30°E 45°E 60°E 75°E 90°E 105°E 

Longitude 


Figure 48. Model’s Sea Surface Elevation Two Hours After the Earthquake 
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Figure 49. Comparison of Jason 1 vs. Model 
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2. Topex/Poseidon 

The second satellite used in this study was Topex/Poseidon. This satellite 
crossed the equator at 3:00 a.m. UTC. Figure 50 presents the model’s sea surface 
elevation two hours and five minutes after the earthquake. Again, the black line 
represents the points used for comparing sea surface elevation in Figure 51 between 
the satellite’s altimetry data (blue line) and the model’s data (red dots). In this 
case, the satellite’s altimetry data contain several gaps. The model estimates very 
accurately the leading wave crest at about 5°S. The gaps in the satellite’s data between 
5°N and the equator do not show the double peak structure that was evident in Jason- 
1 data. Between 12°N to 20°N the model and the satellite’s data are in exceptional 
agreement, and in particular, the model estimates very accurately the leading wave 
crest at about 20°N. 


Sea surface elevation 2 hours 05 minutes after the earthquake 
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Figure 50. Model’s Sea Surface Elevation Two Hours and Five Minutes After the 
Earthquake 
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Sea surface elevation [cm] 


Topex / Poseidon vs. Model 
[2 hrs 05 min after earthquake] 
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Figure 51. Comparison of Topex/Poseidon vs. Model 
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3. Envisat 

Envisat’s altimetry data were also used in this study. Envisat crossed the 
equator at 4:14 a.m. UTC. Eigure 52 presents the model’s sea surface elevation 
three hours and fifteen minutes after the earthquake. The black line, as previously, 
represents the points used to compare the sea surface elevation compared in Figure 
(53) between the satellite’s altimetry data (blue line) and the model’s data (red 
dots). In general, this comparison indicates that the simulation is very accurate in 
representing tsunami waves. In addition to the previous two cases, the reflected wave 
from Sri Lanka between lO'^N to 15°N is reproduced also accurately. The height is 
slightly lower, due to the lack of inundation algorithms, which affects the wave heights 
by the introduction of artificially deep waters at the shallow depth near coastlines. 


Sea surface elevation 3 hours 15 minutes after the earthquake 
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Figure 52. Model’s Sea Surface Elevation Three Hours and Fifteen Minutes After the 
Earthquake 
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Figure 53. Comparison of Envisat vs. Model 
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C. COMPARISON TO FIELD MEASUREMENTS 
1. Mercator 

During the Indian Ocean tsunami, the Mercator, was anchored about 1.6 kilo¬ 
meters off the Phuket coast (Thailand) and specihcally at the location 07.715°iV, 
098.28°i7. Based on the yacht’s depth gauge measurements, the hrst wave (trough) 
arrived at its location at 02:38 a.m. UTC, one hour and thirty-nine minutes after the 
earthquake, while the three main tsunami waves had trough-to-crest wave heights of 
6.6, 2.2, and 5.5 meters. Figure 54 depicts the simulated elevations and also indicates 
the observed arrival time. The model appears to simulate very accurately the arrival 
time of the first wave at the Mercator’s location and that this was negative (wave 
trough), indicating the subsidence on the eastern side of the initialization area. How¬ 
ever, the amplitude is underestimated. The reason for that, by taking into account 
similar difficulties of other models [31], is most likely caused by misrepresentation of 
the water depth. The simulated water depth (grid depth) at the Mercator’s location 
was 45 meters in comparison to an observed of 12 meters. 
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Arrival Time 

"Mercator"s fishfinder vs. Model 



Figure 54. Comparison of arrival times - Mercator’s fishfinder vs. Model 
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D. SYNOPSIS 

From the examination of the foregoing results it is evident that: 

1. The model accurately predicts the positive (wave crest) first wave at all tide 
gauge stations positions due to the uplift of the western side of the initialization area, 
and the negative (wave trough) at the location of the Mercator due to the subsidence 
on the eastern side. 

2. The model simulates very accurately the arrival times and the amplitudes 
of the hrst tsunami waves for the tropical Indian Ocean (Male, Gan, Diego Garcia, 
Port Louis, Salalah, and Pointe La Rue) and the southwest Indian tide gauge stations 
(Okha and Mormugao). All these stations received tsunami waves directly from the 
source area. 

3. The differences in arrival times or in maximum amplitudes between simu¬ 
lation results and records from the Indian tide gauge stations (Paradip, Vishakhap- 
atham, Ghennai, Tuticorin, Kochi, Golombo and Hanimaadhoo (Maldives)) are due 
to the absence of inundation algorithms. According to [22], the tsunami waves that 
arrived at these stations, combined of waves traveling parallel to the earthquake 
rupture and waves reflected from the southwestern coasts of the Indonesian Islands. 
Additionally, in the Hanimaadhoo, Golombo, Tuticorin, and Kochi stations’ cases, the 
tsunami waves had to pass over areas with depths less than 11 meters. Due to the 
boundary conditions used at these simulations (all positions with depths less than 11 
meters were treated as having depths of 11 meters) the model predicts earlier arrival 
times. 

4. The model underestimates the amplitude of the arrived tsunami waves at 
the Mercator’s location due to misrepresentation of the water depth by the used mesh 
grid. 

5. By taking into account that the tsunami wave amplitude squared is propor¬ 
tional to the potential energy, the results depict the west and southwest propagation 
of most of the energy due to the orientation of the earthquake rupture (see Figure 
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55). Only a fraction of the energy propagated toward the south and southeast, which 
is in agreement with [22], 


Maximum heights of waves in meters 



Longitude 

Figure 55. Energy Propagation of the Simulated Tsunami Waves 


6. According to [22] the mid-ocean ridges played a major role as wave guides 
that transferred the tsunami energy to far-held regions outside the source area in the 
Indian Ocean. The model reproduced this feature very accurately as illustrated in 
Figure 55 in combination with Figure 1. Finally, Figure 55 presents, the areas where 
the tsunami focused most of its energy. 
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V. CONCLUSIONS AND 
RECOMMENDATIONS 

This study, as part of a research effort regarding the development of a triangu¬ 
lar discontinuous Galerkin oceanic shallow water model, formatted real bathymetry 
data of the Indian Ocean and simulated the propagation stage of the Indian Ocean 
tsunami that occurred on December 26, 2004. The model’s main advantage is related 
to the geometric (grid) flexibility, which made it possible to represent the tsunami’s 
propagation with different resolutions throughout the Indian Ocean. 

The model’s results were compared to the records of hfteen tide gauge stations 
around the Indian Ocean, satellite altimetry (Jason-1, Topex/Poseidon and Envisat 
data), and held measurements. These results showed that the model is capable of 
producing accurate estimates of arrival times at distant locations. Some considerable 
differences were noticed from the above comparisons due to inhuences by rehections, 
misrepresentation of the water depth by the mesh grid, and the no-hux boundary 
conditions serving as a proxy for inundation algorithms. 

The thesis also derived the analytic solution of the linear Munk problem (two- 
dimensional problem) for future use in validating the diffusion operators that were 
not used in the tsunami simulations. 

Future studies should involve the following: 

1) The use of inundation algorithms in order to complete the simulation of 
the Indian Ocean tsunami December 26, 2004 (propagation and run-up stages), by 
using this DG model. These results could be compared against the already existing 
de-tidal records. 

2) The simulation of the Indian Ocean tsunami December 26, 2004 (prop¬ 
agation and run-up stages), by using this DG model and more complicated initial 
conditions (for example more than two planes). 

3) The simulation of earthquakes with different magnitude for the same region 
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of the Indian Ocean, or the simulation in other parts of the world, in order to identify 
the regions that are likely to be severely affected by tsunami waves, and to establish 
seismic criteria for issuing tsunami warnings in the case of actual tsunamis. 

4) The effects of storm surge modeling (propagation and run-up) caused by 
hurricanes, cyclones and typhoons on coastal areas. 
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APPENDIX A. DERIVATION OF ANALYTIC 
SOLUTION FOR THE MUNK PROBLEM 


In Stommel ’s theory, the inclusion of the bottom friction raises the order of 
the governing equation and enables the beta effect to produce a western boundary 
current [19], which makes Stommel’s solution not quite realistic. The ocean currents 
are concentrated in the upper kilometer of the ocean, they are not barotropic, and 
they are not independent of depth. 

Walter Heinrich Munk noted that the frictional effect due to small-scale eddies 


is crucially important in ocean dynamics, and solved the problem by considering the 
lateral eddy viscosity to explain the ocean current distribution, based on the observed 
wind system. 

Consider a flat-bottomed ocean of depth H driven by wind stress at its surface. 
The momentum shallow water equations in two-dimensions, where bottom drag has 


been neglected while the lateral eddy viscosity is included, are: 

1 dp d'^v d'^v 1 dry 

Pody^ ^ dx^ dy"^ ^ Po dz 

1 dp .d‘^u d'^u, 1 dTx 

Podx^ ^ dx"^ dy"^ po dz 


(A.l) 

(A.2) 


Pressure can be eliminated by cross differentiating these equations and taking the 
difference J|(A.l) — J^(A.2) : 


du dv df 


du. 1 <9 , 

,dv du. Id. , . 


(A.3) 


The use of the continuity equation fy + |]^ + ^ = 0 and beta-plane approximation 
/3 = into Equation (A.3), gives the following: 


/( 


dw 

dz 


+ pv 
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dz 
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dx 
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{curlf) 

{curlf) 


(A.4) 



The Equation (A.4), after multiplying by po, integrating from the bottom (-H) to the 
surface (zero) in the vertical with boundary conditions tc = 0 at 2 : = 0 and = —H, 
assuming zero bottom friction, and using the definition of the streamfuntion: 


gives the following: 


Po 


I tidz — Mx — 
l-H ay 

. < 9 ^ 

/ vdz = My = — 

l-H OX 


= 

Ox 


Coriolis 


z/V^^ + curifo 

LateralFriction Winds 


(A.5) 

(A.6) 

(A.7) 


According to Munk, in the western boundary layer (WBL) Coriolis force is 
being balanced by the lateral friction while at the interior, Coriolis is being balanced 
by the wind’s effect: 


Western Boundary Layer 




Interior 

/ 3 || cuHto 


Consider the idealized model for the North Atlantic subtropical gyre, where: 

- Lx = 5000A;m 

- Ly = 2500/cm 

- Lm is the width of the WBL, {Lm -C L^) 

- Tox = —A COS ^ is the x-component of the wind force 

- 'Toy = 0 is the y-component of the wind force 

- A = 0.1^. 

m 


i. Interior Solution 

The solution of the equation 

= curlTo 
Ox 
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is as follows: 


where 


'h (x, y) = ^ {x) sin 


Try 


T / \ I T \ 

'J' (a^) = ^ - a;) 

Thus, at the interior, < a; < and 0 < y < L^,, the streamfunction solution can 
be written as 


'J'(a^,l/) = ^ -a:)sin (^) (A.8) 

The use of Equation (A.8) and Equations (A.5) and (A.6) of the streamfunction’s 
dehnition gives 


u = 


Att^ 





(A.9) 


V = 


An . (Try 
sin — 


(A. 10) 


[^LyPoH 

The momentum shallow water equations in two-dimensions can also be written as 
follows: 


dh 

,dA 

dA 

-f/yr + 
dy 

'dx'^ 

dy'^ 

dh 

,d^u 

di^u 

ox 

'dx"^ 

dy"^ 


— A cos 


ny 


(A.ll) 
(A, 12) 


where h{x, y) is the free surface elevation. The substitution of the Equations (A.9),(A. 10) 
into (A. 11),(A. 12) gives the following: 


hint{x, y) = {L^ - x) sin 

gfJLyPoH 


cosf^Uc, (A.13) 


/ 


LyJ g(3LlpoH 

where Ci is an unknown constant, and multiplying by the gravitational acceleration 
g, gives the geopotential height solution: 

(Tiy^ 


a) - mxH 


cos 


+ 5'C'i (A-14) 
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Ai = 0 



Since the WBL solution cannot grow exponentially, it must be that C 2 = 0, and 
the three remaining unknowns Ci, C 3 , and C 4 can be evaluated by using the following 
conditions: 


(1) = 0 ^ ^U=o = 0 

(2) = 0 ^ f U=o = 0 ^ f U=o = 0 

(3) at distance Lm the two solutions (WBL and interior) must be the same 
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The use of the above three conditions gives the following: 


T (x) = 


AttL^ 

(3L 


y L 


/ ^/SX 1 VSx' 

1 — e COS —-1-^ sm —— 

V 2Lm Vs 


Thus, at the WBL, 0 < a; < Lm and Q < y < Ly, the streamfunction solution can be 
written as 


T {x,y) = 


AttL^ 


X ( VSx 1 VSx' 
I — e 2 Lm I cos -1- TP sm 


2Lr 


Vs 2Lr. 


sin ^ (A.15) 


As for the interior solution, the use of Equation (A.15) and Equations (A.5) and (A.6) 
of the streamfunction’s dehnition gives the following: 

/ _/ VSx Vs . VSx\\ f7ry\ 


u = 


(A.16) 


2VSAtiL^ _x . VSx . (TTy\ 

V = --- e 2 im sm - sm — A.17 

SpoH /3LyLm 2Lm \^y / 

With the same way as for the interior solution, the substitution of the Equations 
(A.16),(A.17) into (A.11),(A.12) gives 
AnfL^ 


hwbi Vi y') 


gPlyPoH 


. _/ . VSx r- VSx 

1 -e 2 Lm sm-h VO cos- 

3 I 2Lm 2Lm 


AvLx 

g/3poH 


(2VSti^ y3\ . VSx 1 73: 

[siViVn^ sLi) 


X 


. TT-U 

sm + 

Ly 

Try 


e 2im cos -—h C 2 (A. 18) 


where C 2 is an unknown constant, and multiplying by the gravitational acceleration 
g gives the geopotential height solution: 

Att/L^ 


(j^wbi Vi y') 


AuL^ 


LyPtyH 
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